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At non-zero temperature and when a system has low-lying excited electronic states, 
the ground-state Born-Oppenheimer approximation breaks down and the low-lying 
electronic states are involved in any chemical process. In this work, we use a 
temperature-dependent effective potential for the nuclei which can accomodate the 
influence of an arbitrary number of electronic states in a simple way, while at the 
same time producing the correct Boltzmann equibrium distribution for the electronic 
part. With the help of this effective potential, we show that thermally-activated low- 
lying electronic states can have a significant effect in molecular properties for which 
electronic excitations are oftentimes ignored. We study the thermal expansion of 
the Manganese dimer, Mn2, where we find that the average bond length experiences 
a change larger than the present experimental accuracy upon the inclusion of the 
excited states into the picture. We also show that, when these states are taken into 
account, reaction rate constants are modified. In particular, we study the opening of 
the ozone molecule, O3, and show that in this case the rate is modified as much as a 
20% with respect to the ground-state Born-Oppenheimer prediction. 
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I. INTRODUCTION 



In the ab initio study of the behaviour of molecular systems, it is common to perform the 
Born-Oppenheimer separation of the nuclear and electronic parts of the molecular wave- 
function. This approximation is based in the large difference between the masses of the 
electrons and the nuclei,- 1 ^ and therefore becomes exact only in the classical limit for the 
nuclei. Computing the electronic energy spectrum for different positions of the nuclei, one 
obtains the different so-called Potential Energy Surfaces (PESs) — one for the electronic 
ground state, and one for each excited state. Now the question arises, how may we use 
these PESs in order to produce accurate and convenient physical models? The simplest op- 
tion is ground-state Born-Oppenheimer (gsBO), or typically just Born-Oppenheimer (BO). 
One may then consider the nuclei to be quantum particles and solve their corresponding 
Schrodinger equation, or take the classical limit and perform BO Molecular Dynamics. In 
any case, a model based on only one PES (usually the ground state one) is an adiabatic 
approximation, based on the neglect of the non-adiabatic couplings. 

However, at non-zero temperature and when a system has low-lying excited electronic 
states, these electronic states are involved in any chemical process, and their influence pro- 
duces the so-called non-adiabatic effects. In this paper, we use a thermodynamically ac- 
curate generalization, introduced in Ref. [j, of the gsBO potential, built with the PESs of 
the lowest-lying electronic states appropriately weighted to produce the correct Boltzmann 
equilibrium distribution of the electronic part, to specifically study the thermal expansion 
and reaction rates of small molecules. For both phenomena, a substantial amount or re- 
search is constructed on the implicit assumption of one electronic potential, which may be 
fitted to experimental results or computed ab initio. The possibility of electronic excitations 
is typically either ignored or handled by independently considering each of the excited PES 
(although we can mention Ref. |4| as a very recent exception to this general trend, with simi- 
lar aims to the ones we pursue here). The thermally averaged potential energy surface that 
we use in this work permits to include electronic excitations while still preserving the single 
potential methodology. 

The possibility, in which our formalism is based, of dividing a system of particles into 
a quantum and a classical subsystem (typically, electrons and nuclei) is of wide interest in 
several areas of physics and chemistry. If the temperature is of the order of the electronic gap 



or larger and excited electronic energy levels have to be included in the formalism, a variety 
of approaches can be considered.- - — For example, the decay-of-mixing method by Truhlar 
and coworkers^ constitutes a powerful mixed quantum-classical scheme for modeling non- 
Born- Oppenheimer chemistry, although the incorporation of temperature to these methods 
has not been studied as far as we are aware. In the case of Ehrenfest dynamics, which 
also includes non-adiabatic effects at a different level, the temperature has been introduced 
through the formulation of Nose.— ^ 

The temperature dependence of molecular properties (geometry and vibrational frequen- 
cies) of free molecules has been a subject of research for more than 60 years,— with new 
theoretical analyses coming out still in very recent times.— ^ From the experimental view- 
point, on the other hand, several studies of hot molecules^— have found thermal expan- 
sion of bond lengths. In this paper we calculate the temperature dependence of the bond 
length in the case of diatomic molecules in which nuclei can be treated as classical par- 
ticles. Although bond-length expansion as the temperature increases is expected even for 
temperature-independent potentials such as the gsBO one, the use of our thermodynam- 
ically more accurate effective potential adds new non-adiabatic effects which, as we will 
show for the Mn 2 dimer, may significantly alter the final quantitative results if low-lying 
excited states are present. These effects must be considered when performing the necessary 
rovibrational averaging in order to compare accurate theoretical and experimental results. 

The temperature dependence of the transition rate is also a constant subject of study22r— 
and, more recently, the question has been asked of whether or not quantum tunneling below 
the energy barrier significantly enhances the reaction rate.— ~— Using a general framework 
to describe tunneling,— it is shown that tunneling below the barrier only occurs for temper- 
atures less than a reference one, denoted by To, and which is determined by the curvature 
of the ground-state PES (gsPES) at the top of the barrier. However, at non-zero tempera- 
tures and when a system has low-lying excited electronic states, all estimates based on the 
ground-state PES should be reconsidered. 

As demonstrated in this work using our thermodynamically accurate generalization of the 
gsBO potential, the inclusion of low-lying electronic states into the picture may significantly 
alter the reaction rates and the curvature near the top of the barrier. As a model system to 
exemplify the added effects, we have chosen the transition between the open and closed forms 
of the ozone molecule, in which the barrier lies in a region where there is an avoided crossing 
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between the ground and first-excited electronic states. Having no problems whatsoever with 
avoided-crossing situations, our effective potential is a convenient choice to account for the 
influence of low-lying excited states in this phenomenon at non-zero temperature. 

In Sec. [TT1 we provide a comprehensive summary of the definition and meaning of the 
effective potential introduced in Ref. O which includes the effect of excited electronic states 
and which will be used throughout the document. In Sec. IHI A\ we present the first ap- 
plication of the effective-potential technique: the study of the thermal expansion of the 
Manganese dimer, Mn 2 ; in Sec. IIIIB[ we discuss the influence of the inclusion of low-lying 
excited electronic states on the reaction rate of small molecules. Reaction rates are affected 
on the one hand by the energy barrier reduction. On the other hand, the curvature at the 
top of the barrier is smaller for our effective potential than for the gsPES. In this section 
we will also study the case of the opening of the ozone molecule, O3. Finally, in Sec. IIVI 
we comment on the most important conclusions of this work and highlight some possible 
implications and future lines of research. 

II. THEORETICAL FRAMEWORK 

Let us have a quantum-classical system formed by N classical particles described by their 
Euclidean coordinates R := (R±, . . . , Rn) and momenta P := (Pi, . . . , Pn), and n quantum 
particles described by an n-body wavefunction \ip). The starting point of our formalism is 
the assumption that the following is an accurate formula to compute canonical equilibrium 
expected values of quantum-classical observables 0(R,P): 



where k-Q is the Boltzmann constant and T the temperature, dRdP denotes integration over 




(2.1) 



all position and momenta in the appropriate ranges. The object H (R, P) is the quantum- 
classical Hamiltonian, which, in the case of molecular systems, takes the form 




(2.2) 



where 1 denotes the identity matrix, is the mass of the K-ih nucleus, and the electronic 
Hamiltonian, H e (R), contains all particle interactions and the electronic kinetic term (see 
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Ref . [2I for an explicit expression) . 

It is also convenient to write Eq. (12. D as 



(6(R,P)) = J dRdP\x(6{R,P)p^(R,P)) , (2.3) 
in terms of a (R, P)-dependent equilibrium density matrix, defined by 



H(fl,P) 



MR, p) ■= 7 ■ ^ 

dR'dP' tr e k B T 



As shown in Ref. [j, the justification of Eq. (12. ip for computing equilibrium expected values 
in quantum-classical models stems on the one hand from plausibility arguments related to 
the classical limit of the behaviour of the nuclei, such as the ones given in Ref. |27|. However, 



it can also be obtained as the zero-th order approximation (in the quantum-classical mass 
ratio) to the canonical equilibrium associated with an, in principle, more rigorous quantum- 
classical formulation based on the Wigner formalism, as shown by Kapral and Ciccotti^ and 
Nielsen et al.— . It may also be rationalized by entropic arguments^. In any case, and as it 
can be seen in the references in Ref. [3], irrespective of how good an approximation Eq. (12. ip 
is for any given application — always a difficult question — , it is often the desired, target 
expectation value when designing quantum-classical schemes. 

The main realization in which the effective potential that we will use in this work is based 
is that, for observables which do not depend explicitly on the electronic degrees of freedom, 
i.e., which are of the form 

6(R,P) = 10(R,P), (2.5) 
where 0(R, P) is a number, the target expected value in Eq. (12. ip can be rewritten as 

I dRdPO(R,P)e k B T 

(0{R,P))= J - ? • M 

dRdPe k B T 



Now, H eS (R, P; T) is a purely classical, T-dependent, effective Hamiltonian defined as 
H cS (R,P;T) := -k B T lntr e'^^ 

N fin 

EPl _Heim 
— — -k B Thitie k * T 
2M K 

K=l A 



where we have used Eq. (12.21) . In the last line, we have finally defined the promised, purely 
classical, T-dependent, P-independent, effective potential 

He(fl) 

V cS (R;T) :=-k B T Intie (2.8) 

which can be used to describe the behaviour of the nuclei in a classical mechanical set- 
ting, producing the correct target equilibrium in Eq. (12.11) for classical observables — by 
construction — , and including the influence of all the electronic excited states. 

Indeed, if we consider the adiabatic basis {\tp k (R))}, which diagonalizes the electronic 
Hamiltonian H e (R) for each fixed position of the nuclei, 

H e (R)\MR)) = E k (R)\ip k (R)} , (2.9) 

being {E k (R)} the corresponding PESs, i.e., the eigenvalues of the electronic Hamiltonian 
as a function of the nuclear positions, then we can rewrite the effective potential in Eq. (12.81) 

as 

V eS (R;T) = -AfeTln^e k * T (2-10) 



= E (R)-k B T\n 

where we have defined 



e 



k-aT 



k>0 



AE k0 (R) := E k (R) — E (R) . (2.11) 

The expression in Eq. (I2.10p explicitly shows the difference between the ground state 
PES, E (R), i.e., the gsBO potential, and the effective potential V e s(R;T) introduced in 
Ref. y and used in this work. In particular, it is worth remarking that 

• At T = 0, our effective potential becomes the gsBO one. Indeed, it is easy to see from 
Eq. (12TT0|) that 

Urn V cG (R;T) = E (R) , \/R. (2.12) 

• The same reasons that produce the previous result allow us to include in the sum in 
Eq. (I2.10p only the lowest-lying electronic states and still get a good enough approxi- 
mation to the exact expression if the temperature is low compared to the states we are 



neglecting, i.e., if k-^T <C AE^i^R) for them. This fact will be used in the practical 
cases presented in the next sections. 



It can be seen from Eq. (I2.10P that an exact property of the effective potential is that 
it is strictly lower than the gsPES, i.e., V e g(R;T) < E (R), Vi?, T. However, since an 
additive constant in a potential energy is not measurable, it must be noticed that this 
inequality is relevant only inasmuch the difference Eq(R) — V e s(R; T) actually depends 
on R. See Sees. IIII Al and IIIIBI for concrete examples of this situation. 

As we discussed in Ref. [j, from Eq. ( I2.1(jp . we see that, if the second derivative of 
AEiq at a barrier top qs verifies 



8 2 AE 
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dq 2 



(Qi 



> 



1 / dAE 



knT 



(Qi 



then we can prove 



d 2 V r 



off 



dq 2 



< 



d 2 E 



dq 2 



(2.13) 



(2.14) 



i.e., the curvature of the effective potential at the barrier top is smaller than the one 
associated to the gsPES. In avoided crossings — a very interesting general case, and 
the one actually studied in Sec. IIIIBI — since the first excited state approaches the 
gsPES and then recedes from it, we have that the derivative (OAEiq/ dq)(qB) will be 
approximately zero and the condition in Eq. (12.131) will be approximately satisfied, 
together with Eq. (12TT4")) . 

Note that in order to obtain the effective potential, and the corresponding averages, it is 
not necessary to compute the non-adiabatic couplings. This is a consequence of the fact that 
the purpose of the effective potential is the computation of averages at canonical thermal 
equilibrium, and not the dynamics that may lead to it. The non-adiabatic couplings are 
necessary to carry the system to equilibrium - by providing the necessary channels to couple 
the various electronic states. Because of this, we are considering that the difference between 
these averages and the ones that would be obtained by considering the gsPES only is a non- 
adiabatic effect. However, the equilibrium averages predicted by Eq. (12.11) can be obtained 
without explicit consideration of the couplings. The magnitude of those couplings might be 
very relevant to compute the speed of the thermalization: small couplings may necessitate 
long thermalization times, but those analysis are beyond the scope of the present work. 
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In this respect, also note that all excited electronic states are included in the definition of 
the effective potential - because all states are included in Boltzmann's equilibrium formula, 
regardless of how they may be coupled by external fields or, in the quantum-classical case, 
by the non-adiabatic couplings. The weight of each state is entirely determined by its energy. 
Any state is present in the averaging, even if its couplings to the ground state and to any 
other accessible state is zero. This is an example of ergodic difficulties, and obviously, a 
dynamical averaging would not include such a state unless it is already included in the 
initial state sampling. 

The straightforward application of our scheme would then be inadequate if one is inter- 
ested in computing a "restricted equilibrium" average, in which a state (or full subspace of 
states) is known to be absent, due to symmetry rules. The "experimental average" would 
not contain those states, even if the true canonical ensemble average does. However in this 
situation it would be easy to correct our scheme by simply not including the forbidden states 
in the formulas. 

III. RESULTS 

A. Non-adiabatic effects on the thermal expansion of the Mn 2 dimer 

Theoreticians usually identify the "molecular structure" with the equilibrium structure, 
i.e., the point determined by the absolute minimum of the ground state Born-Oppenheimer 
potential energy surface (gsPES). This point in ~R 3N space (A being the number of atoms) 
is well defined and has an easy intuitive meaning: the position occupied by the nuclei at 
equilibrium, if they had infinite mass (in which case they would be classical point particles). 
This geometry corresponds to a motionless molecule, that does not exist because molecules 
vibrate and move even at zero temperature. Therefore, this equilibrium structure is a theo- 
retical concept that is not provided — at least not directly — by the experimental techniques 
utilized to probe molecular structure. 

In fact, different experimental techniques yield different averaged results, whose value 
and meaning depend on the physical process involved in the measurement. For example, 
X-ray diffraction provides distances between the electronic charge distribution centroids. 
Gas-phase electron diffraction (ED) , on the other hand, provides internuclear distances. Mi- 
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crowave spectroscopy measures moments of inertia, which can be directly related to nuclear 
distances to the molecular center of mass. Infrared, Raman, and ultraviolet spectroscopies 
can also be used for complementary analysis. In all cases, the results are averaged over the 
populated rotational and vibrational molecular states — and, as we will show below, possibly 
over different electronic potential energy surfaces. Those techniques achieve a remarkable 
precision (of the order of 0.001 A), but nevertheless provide different numbers each one of 
them. 

In order to compare the results obtained in different experiments and in precise ab initio 
calculations, it is necessary to use a "common denominator" representation, which can very 
well be the equilibrium structure mentioned above, usually called the r e structure. One must 
therefore know how to relate the experimental result to this concept. In an ED experiment 
at a given temperature, for example, one obtains the so-called r a structure, an operational 
concept with no clear physical meaning. It can be related, however, to the thermally averaged 
internuclear distance, or r g structure, which is of physical significance. It is not equivalent 
to the r e structure, not even at K, because of the vibrational and rotational (also called 
centrifugal) distortions. The relationship between r e and an averaged structure such as 
r g is not straightforward, but must be considered if we want to validate high precision 
theoretical ab initio calculations with experimental results or vice versa. This relationship 




FIG. 1. Six lowest lying PESs of the Mn2 dimer, taken from Ref. [311 ] . The calculations were 
performed through multireference variational calculations coupled with augmented quadruple cor- 
relation consistent basis sets. These PESs correspond to the ground state "manifold", that correlates 
to ground state Mn atoms. 
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was first considered by Bartell,— and has later been developed by several authors.— ~— It 
was soon realized that, in general, molecules expand as temperature increases, due to both 
the anharmonicity of the vibrations and to the centrifugal "force" that rotations exert on 
the structure. 

In all these studies, an assumption is implicitly made: the electrons are adiabatically tied 
to their ground state, and therefore the analysis is performed by considering the gsPES only. 
However, as discussed in the introduction, the existence of non-zero non-adiabatic couplings 
permits the system to visit electronic excited states, and the thermodynamic averaging 
should account for this possibility if the energy gap between the gsPES and the excited ones 
is not large in comparison with the thermal energy, k#T. If calculations and experiments 
are to be succesfully compared, it is therefore necessary to consider the possible effect of the 
excited electronic states. We propose the use of our thermally averaged potential energy 
surface introduced in the previous section for this purpose. 

In principle, the rotational-vibrational averaging necessary to compute r g must be per- 
formed assuming quantum-mechanical nuclei. This is obvious at K, where the zero-point 
vibration would be completely absent in a classical treatment. However, we are mostly inter- 
ested in the high-temperature regime — where the excited electronic states are expected to 
play a larger role and the system becomes classical. This fact can be exemplified by looking 
at closed — albeit approximate — theoretical formulae that exist for the simplest cases on 
which we concentrate in this study: diatomic molecules. By truncating the Taylor expansion 
of the PES around the equilibrium bond length, i.e.: 



Toyama et al.— computed, in the quantum case, the temperature-induced variation of the 
average internuclear distance (with respect to r e ) as: 



where u e is the vibrational frequency at equilibrium (yj e = \Jkij '//, where fi is the reduced 
mass of the two nuclei). 

The first term in Eq. (13.21) is the centrifugal distortion; it arises from the harmonic 
potential term in Eq. (13. ip . and can be considered as the variation in the average bond 




(3.1) 



(Ar)(T) := (r)(T)-r e 
k 2 r e 2k\ 



coth 



2k B T 



(3.2) 
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length caused by the overall molecular rotations. Interestingly, a classical treatment leads 
to the same expression for this term. The second term, on the contrary, has a genuinely 
quantum behaviour at low temperatures. In fact, it does not approach zero as T — > K, 
producing a zero-point vibration variation of the bond length with respect to r e . We can 
now take the classical limit by considering fi — > oo, in which case this second term becomes 
linear in T: Sh^k^T/k^. But note that the same behaviour occurs if we take k^T large with 
respect to uj e — i.e., the system becomes classical for large enough temperatures. 

In view of this, and since we are interested in the high-temperature regime in which 
our effective potential may significantly differ from the gsPES, we have used the classical 
approximation to compute the average bond length. To this end, we used Eq. (12. 6p . which 
for the dimer case (and considering that the function 0(R,P) is in this case nothing else 
than the internuclear distance r = \Ri — R 2 \) , reads: 



Note the presence of a somewhat arbitrary upper limit of integration L. This value cannot 
be made arbitrarily large, since in the limit L — > oo, the equilibrium bond length is also 
infinite (assuming that, as it is always the case, the potential function remains finite at large 
internuclear distances). In fact, at equilibrium and at any non-zero temperature, a dilute 
gas of diatomic molecules in an infinite space does not really contain dimers but isolated 
atoms. In the real world, dimers exist because there is always some form of container, or 
they are in a very long-lived metastable state. In practice, one must choose a value of L 
such that the results are almost insensitive to small variations of it — acknowledging that, 
if L is increased to very large values, the value of (r)(T) will start growing. 

The question that we want to answer here is whether or not the use of the effective po- 
tential in Eq. (13.31) leads to significantly different results with respect to the results obtained 
using only the gsPES, i.e.: 



The answer cannot of course be universal, and depends on the chosen system and the 
temperature regime observed. In order to illustrate the issue, we have concentrated on the 
Mn 2 molecule; a Van der Waals weakly-bound molecule and a specially difficult theoretical 
case,— for which a good feedback between experiments and theory could help validate 
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FIG. 2. Calculated variation of the average internuclear distance as a function of temperature 
for Mn2. (Ar)eff is the result computed using, in Eq. 13. 31 the effective potential including all six 
lowest-lying states in Fig.[H while (Ar)o is only computed with the gsPES. (Ar)g is the third-order 
approximation to the quantum result, as calculated using Eq. ()3.2|) . and (Ar)g lass is its classical 
limit. (Ar) arm is the centrifugal term, common to all cases. 

the conflicting theoretical results. As we shall show, the existence of very low-lying excited 
electronic states in the vicinity of the equilibrium distance has an important effect on (r)(T), 
and therefore it is crucial to consider them to make proper comparisons. 

To build our averaged potential denned in Eq. 12. 10} we depart from the potential energy 
curves provided by Tzeli et al.— , computed from first principles with very accurate multiref- 
erence variational calculations. We only consider the six almost-degenerate, lowest-lying 
states, displayed in Fig. [TJ We have adjusted these curves to Morse functions, i.e.: 



By finding the best match for the parameters D,a,r e and Voo, an almost perfect fit can be 
obtained with respect to the results of Tzeli et al. The plot clearly shows the reasons for 
choosing Mn 2 in this study: the lowest-lying states are extremely close to each other, and 
the potential well is very shallow. 

The results are shown in Fig. [2j The two key curves are the ones denoted by (Ar) e g and 
(Ar)o, which are the results obtained with Eqs. (13. 3p and (13. 4p . respectively. The difference 
due to the use of the effective potential instead of the ground state one is notable. It is larger 
than the resolution of modern experimental techniques, even in the lower-temperature range. 
One may then conclude that any assessment of the quality of a theoretical model based on a 



V(r) = D e 



-2e 



+ V C 



oo ■ 
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FIG. 3. Variation of the average internuclear distance, at 300 K, for a fictitious dimer with 
two relevant PES, as schematically shown in the inset. The results are given as a function of the 
displacement of the excited PES with respect to the gsPES, in two directions: variation of the 
energy minimum A£io, and difference in the position of this minimum, Ar e . 

comparison to experimental results should consider the influence of these low-lying electronic 
excited states. 

In Fig. [21 we also display the approximate quantum result given by Eq. (13. 2p , and denoted 
by (Ar)o- This quantum curve is only valid at low temperatures, since it stems from a 
perturbative truncation of the potential. We display it in order to demonstrate how the 
system behaves almost classically in most of the temperature range of the plot, thus justifying 
our classical treatment. Indeed, if we plot the classical limit (/x — > oo) of Eq. (13. 2p . denoted 
by (Ar)Q lass , it quickly becomes almost identical to (Ar)g. In this temperature region, our 
classical calculation, which does not truncate the potential curve, should be almost exact. 
Finally, for completeness, the curve (Ar) arm is the centrifugal term — the difference with the 
rest of the curves would be the vibrational contribution in each case. 

Beyond this particular example, a more general question has to be addressed when must 
one expect the electronic excited states to influence the thermally averaged internuclear 
distances — and therefore any experimental measurement of molecular geometry? A simple 
visual inspection of a few lowest lying excited PES, and a simple calculation with our ther- 
mally averaged PES should give us a quick answer. Two key parameters should be carefully 
examined: the "gap", or difference between the gsPES and the closest excited ones, and how 
much the position of the minima of the two curves differ. This fact is illustrated in Fig. [21 
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where we have considered a fictitious dimer with two closely lying PES. The parameters of 
the gsPES correspond to the Hydrogen molecule, whereas the first excited PES is a rigid 
displacement in two directions: varying the minimum energy (Ai?io), and the position of 
the corresponding minimum (Ar e ). 

The 2D plot displays in Fig. [3] (Ar) c g at room temperature (300 K). As the gap becomes 
large, the results converge towards (Ar) , the thermal expansion entirely due to the gsPES. 
The plot shows how, for the results to differ significantly, the gap should not be larger than 
a few mHa — which is easily understood since, at 300 K, k^T is approximately 1 mHa. But 
also note that, even if the gap is small, there is no change with respect to (Ar)o if the 
positions of the minimum of the two curves do not differ (Ar e is small). In other words, if 
the two PES are merely a rigid vertical shift of each other, the thermally averaged PES is 
also a rigid vertical shift and nothing change. 



B. Non-adiabatic effects on the reaction rate and the opening of ozone 

Reaction-rate theory focuses on studying the behavior at long times of systems with 
different equilibrium states. This is a subject of great interest in many biological, chemical 
and physical problems. As noted by Arrhenius in 1889,— the cornerstone of the theory is the 
e -A/T temperature dependence of the reaction rates. Such dependence can be understood 
in the framework of a transition-state theory where the system evolves as a function of a 
given reaction coordinate q from the metastable state A to the metastable state C through an 
energy barrier B, being the activation constant related with the barrier energy U = Eb — Ea 
(see Fig. gj). 

In general the reaction rate can be written as 

r = ke -U/k B T^ ( 3 g) 

where the prefactor k depends on the temperature T, the friction coefficient or 'damping' 
of the system, and the details of the potential energy function. An accurate estimation 
of this prefactor has been shown to be a formidable task and many articles have been 
devoted to this end — deserving an special mention the celebrated one by Kramers^. See 



also Refs 



23 



43 



and 



44 



As it has been shown in Ref. [3] and summarized in the previous sections, the inclusion of 
excited electronic levels becomes important at certain temperatures for obtaining suitable 
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PESs of molecular systems, different from the gsBO one. The object of this section is to 
consider the effect of these excited states on the thermal-activation rate calculations. 

Let us denote by r e fj the activation rate for the effective potential V e s(q; T) in Eq. (I2.10p . 
and by tq the rate for the gsPES, Eo(q) (both of them expressed as a function of the 
reaction coordinate q). The main effect of the inclusion of the new terms associated to the 
excited electronic states is a reduction of the energy barrier from Uq = E^qs) — E (q A ) to 
U c s = V eS (qf) - V cS (qf). We will have in general qf ~ q A and V eS (qf) ~ E (q A ), and 
thus the change in the energy barrier is given by the change at the potential maximum. 
Regarding the activation rate, in the simplest picture, we find that r e fr ~ e -u c e/k B T an j 
fcfr/?"o ~ e ~ (Uch -Uo) / k B T ^ Therefore, if U c q — Uq is large enough, the effect of the energy 
barrier reduction on the activation rate will be important. In addition to this effect, it 
is also important to realize than the activation rate will also show a deviation from the 
expected e~ u ^ kBT temperature dependence due to the fact that U e ff itself depends on T. 
With this in mind, the importance of the excited electronic states can be studied looking 
for deviations of r(T) from its expected dependence. 

A more detailed analysis must take into account the influence of the prefactor k too. To 
this end, some approximations need to be done. Let us place our problem in the context of 
the large-barrier (U/k&T 3> I) and strong-friction limit. In such a situation, we learnt from 




q A q B qc 



FIG. 4. Sketch for the usual reaction rate problem, in terms of a generic reaction coordinate q: Two 
metastable states at q A and qc are connected by a barrier with height U = Eb — E A and maximum 
at qB- The curvatures are also shown. 
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FIG. 5. Calculated potential energy profiles for ozone. Top: Eq((/)), Ei{4>) and V c s(4>) at 500 K. 
State A lies at (pA = 60° {Eq^a) = K) and the barrier is located at approximately <pB = 83.4° 
{Eq{(Pb) = 9094 K). Bottom: Potential profiles close to the barrier. Eo((p), E\{(j)) and V e s(4>) for 
T =100, 200, 300, 400 and 500 K. 



Kramers that,— 



^khd = — 

7 



2vr 



JJ(q)/k B T 



dq 



' 1A 



(3.7) 



where U(q) := V(q) — V(qA), being V(q) the appropriate potential. Here, uja is associated 
to the curvature of the potential at point A, 7 is the friction coefficient and KHD stands for 
Kramers high-damping limit. 

This formula can be used directly to calculate the reaction rate by performing the integral 
numerically, and this is what we will do in this section. However, before doing that, let us 
introduce two simple approximations that are instructive and give some insights about the 
general features of the problem. For a large barrier, a narrow region around the maximum 
gives most of the contribution to the integral. In many cases we can use the so-called 
parabolic barrier approximation^: U ne3iVB ~ Ub — | w s(<? — Qb) 2 - Then 



KHD 



2?T7 



In this situation, we have 



Tcff ^ uf(T) e -(u eS ( T )-u )/k B T 
r ub 



(3.8) 



(3.9) 



Now, if Kg ~ Eo beyond the barrier and V e s < Eq at the barrier, we expect uj c g < ujb 
(see also the end of Sec. IIIip . Hence, the prefactor effect opposes the exponential one: The 
rate will typically diminish because of the changes in the prefactor, but it will typically 
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FIG. 6. r c s/ro ratio as a function of T (red circles) as computed numerically using Eq. ()3.7[) . We 
also show the value of the exponential reduction factor e s ( T ^ ksT (blue squares) and the prefactor 
change b(T) estimated from Eq. (|3.13p (green triangles). 

increase due to the changes in the barrier. In any case, given its exponential dependence, 
the effect of the barrier reduction is expected to be dominant in most cases. 

Special care must be taken if the region of the potential close to the barrier (i.e., the one 
that contributes significantly to the integral in Eq. (I3.7p ) cannot be accurately approximated 
by a parabolic function around its maximum. Another common approximation to compare 
with is given by a cusp barrier: 

{U B + mi(q-q B ) if q < q B 
(3.10) 
U B -m 2 (q-q B ) if q > qs 

where m 1; m 2 > are the slopes at each side of the barrier. In this case, the activation rate 
can be written as 

^%rM Xe ~ w ' (311) 

where A := m\m 2 j (mi + m 2 ), and we have 

r e $ _ A e (T) c _ (u ^ (T) _ Uo)/kBT ^ , 312 , 
r A 

In conclusion, we can write in both cases 

T -^~b(T)e 8 ^ k * T , (3.13) 
r o 

where S(T) := Uq — U c g(T) > is the barrier reduction, and b(T) accounts for the changes 
in the prefactor. For the simplest case where E (q) and V e g(q) show a maximum at the 
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same coordinate q B (the position of the maximum is not affected by the new terms), it is 
easy to see that e 5 ^ T ^ kBT ~ (1 + e" A ^ BT ) with A := Ei(q B ) — E (q B ), where we have only 
considered the first excited state E\. 

As a working case, we will consider now the case of the ring opening reaction of ozone, 
which has been previously studied in Ref. ^. Fig. [5] shows the potential profile for ozone as 
a function of the opening angle <fi at different temperatures. Since this is enough for our 
purposes, we consider only the gsBO PES, E (<f)), and the first excited state one, Ei((j)) 
— taken from Ref. [j, where they where computed using the CASSCF method (complete 
active-space self-consistent field).— The effective potential V e s(4>) is constructed with them. 
As we can clearly see in the bottom graph, the potential profile is modified by the inclusion 
of the new term corresponding to E\ (0) only in a narrow region close to the barrier. Also, in 
this case, the potential barrier is close to 9000 K and thus U/k^T ^> 1 for the temperature 
range of interest, which justifies the basic approximations behind the formulae in this section. 
However, there exist a clear asymmetry between E (<f)) and Ei(<f)). As a consequence, the 
barrier location moves with T, and the barrier energy reduction S(T) does not follow the 
expected T dependence. In addition, as we can see in Fig. [51 neither the parabolic nor the 
linear cusp approximations will be suitable to fit the barrier profile close to the maximum. 
Hence, we have directly used Eq. (13. 7p to numerically compute the r e g/r ratio of the system. 
The results are presented in Fig. where we can see that, in spite that the small barrier 
reduction observed, the activation rate changes up to a 20% upon the inclusion of the 
excited electronic states. This indicates that the non-adiabatic effects associated with low- 
lying states should be included in any analysis of this kind if one aims for high accuracy 
in the predictions. Besides, in the same Fig. [6J we also plot the contributions to the ratio 
r cs/ r o coming from the changes both in the potential barrier heigth and shape. As expected 
the increase of the rate motivated by the barrier reduction is moderated by the prefactor 
change. The two effects work against each other, and the exponential dependence on the 
barrier overcomes the influence of the curvature. 



IV. CONCLUSIONS 



In this work, we have shown that the ground-state PES, as defined in the Born- 
Oppenheimer approximation and typically used for many applications in chemistry, physics 
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and materials science, is not the only electronic state that significantly contributes to the 
theoretical prediction and calculation of thermodynamic observables at non-zero temper- 
ature already for small molecular systems. Although this fact was probably expected by 
the reader, we provide actual numerical examples of relevant observables being significantly 
modified by the inclusion of thermally-activated low-lying excited electronic states in real 
systems at not-too-high temperatures: The average bond length of the Manganase dimer is 
shown to change in an amount which is accessible to modern experimental techniques, while 
the reaction rate of the ring opening of ozone is shown to change up to a 20%. Moreover, our 
compact effective potential, which includes any number of such states and which, despite 
its simplicity, is able to produce the correct Boltzmann weight for the electronic subsystem 
— a long sought property in the field of quantum-classical models. 

Also, as discussed in section Sec. Ull a general result when using our effective potential is 
that, in any avoided crossing situation, the curvature on the top of the barrier is smaller than 



the gsPES curvature. As shown in Ref. 



25 



if the curvature is smaller, the tunneling below 
the energy barrier will occur at lower temperature. Therefore, the inclusion of low-lying 
electronic states is important if one wants to adequately discriminate possible quantum-like 
behaviour of the nuclei from simple classical effects due to the direct influence of temperature 
on the potential — for example in biological systems.— 1 ^ 

Additionally, although the effective potential used in this work has been derived in Ref. O 
assuming classical nuclei and equilibrated electrons, it could also be used as an external 
potential to perform calculations on quantum nuclei. This procedure would allow to study 
low temperatures, while still including a possible correction due to electronic excitations. 
Notice, however, that, although the use of the gsPES in the BO approximation as a potential 



for quantum nuclei can be rigorously justified (see, e.g., Eq. (17) in Ref. 1461 ). in the case of 
the effective potential used in this work, its use in this manner is not justified in principle 
because of its intrinsic non-adiabatic origin. 

Finally, it is also reasonable to expect that the use of our temperature dependent effective 
potential provides new insights that could lead to answer the intriguing question of the 
thermodynamical stability of some diatomic dications,— an issue we plan to address in the 
next future. 
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